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ABSTRACT 

We have computed models of rotating relativistic stars with a toroidal magnetic field and in- 
vestigated the combined effects of magnetic field and rotation on the apparent shape (i.e. the 
surface deformation), which could be relevant for the electromagnetic emission, and on the 
internal matter distribution (i.e. the quadrupole distortion), which could be relevant for the 
emission of gravitational waves. Using a sample of eight different cold nuclear-physics equa- 
tions of state, we have computed models of maximum field strength, as well as the distortion 
coefficients for the surface and the quadrupolar deformations. Surprisingly, we find that non- 
rotating models admit arbitrary levels of magnetisation, accompanied by a growth of size 
and quadrupole distortion to which we could not find a limit. Rotating models, on the other 
hand, are subject to a mass-shedding limit at frequencies well below the corresponding ones 
for unmagnetised stars. Overall, the space of solutions can be split into three distinct classes 
for which the surface deformation and the quadrupole distortion are either: prolate and pro- 
late, oblate and prolate, or oblate and oblate, respectively. We also derive a simple formula 
expressing the relativistic distortion coefficients and that allows one to compute the surface 
deformation and the quadrupole distortion up to significant levels of rotation and magnetisa- 
tion, essentially covering all known magnetars. Such formula replaces Newtonian equivalent 
expressions that overestimate the quadrupole distortion by about a factor of five and are inad- 
equate for strongly-relativistic objects like neutron stars. 
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1 INTRODUCTION 

In addition to rapid rotation, also strong magnetic fields can intro- 
duce significant deformations in neutron stars, as shown, for in- 
stance, by |Bocquet et al.j (1995) and Cardall et al. (2001 1, who have 
computed fully nonlinear models of relativistic stars with a poloidal 
magnetic field. At the end of the collapse of the core of a mas- 
sive star, differential rotation could create strong toroidal magnetic 
fields of the order of 10 16 G - 10 17 G inside the hot proto-neutron 
star ( |Bonazzola & Haensel|1993HBonanno et al.|2003||Naso et al.| 
2008 ). As a result, realistic models of magnetised relativistic stars 
require the simultaneous presence of both poloidal and toroidal 
field components. As pointed out in Gourgoulh on & Bonazzolaj 
{1993}, however, this requires a formalism capable of dealing with 
non-circular spacetimes (i.e. with convective currents in the merid- 
ional planes) like the one presented here, but which has so far 
been implemented only in a perturbative scheme ll oka & Sasaki| 
|2003| [2004 >. Triggered also by the increasing interest in strongly- 
magnetised neutron stars due to their relation with soft-gamma re- 
peaters and anomalous X-ray pulsars ( |Duncan & Thompson|1992) 
Thompson & Duncan 19961, a growing number of studies adopt- 
ing perturbative techniques have appeared investigating either the 
field geometry and neglecting the influence of the magnetic field on 



the matter distribution (Ciolfi et al. 2009), or solving the coupled 
Einstein-Maxwell-Euler system, from which the magnetic defor- 
mation can be calculated ( Ioka & Sasaki 2004 ; Colaiuda et al. 2008 ; 
ICiolfi et al.|2TTiO[|Gualtieri et al.|201 l||Yoshida et al.|2012) . Until 
recently, however, fully nonlinear models of relativistic magnetised 
stars were restricted to purely poloidal magnetic fields ( |Bocquet| 
et al. 1995 ; Cardall et al. 20011, for which the generated spacetime 
is circular like in the unmagnetised case (Carter 19731. Following 
the recent insight that also a magnetic field with only a toroidal 
component is compatible with the circularity of spacetime ( |Oron| 
|2002[ >, studies of relativistic models of rotating stars with a toroidal 
magnetic field h ave emerged jFrieben & "R ezzolla 2007; Kiuchi 
|& Yoshida|2008l|Kiuchi et al.|2009||Yasutake et al.|2010||27)TTT ^ 
These new studies have complemented earlier Newtonian investi- 
gations lSinha|1968||Sood & Trehan|1972[|Miketinac|1973^ , which 
being simpler, allowed for the investigation of more complex field 
geometries, in particular of mixed poloidal and toroidal magnetic 
fields (Tomimura & Eriguchi 2005 , Yoshida et al. 2006 ; Haskell 
|et al.|2008[|Lander & Jones|2009||Fujisawa et al.|2012| >. 

In this work, we are mainly concerned with the deformation 
of neutron stars endowed with a toroidal magnetic field because 
they might be an important source of gravitational radiation due to 
the prolate deformation induced by the magnetic field and provided 
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that the axes of symmetry and of rotation are different (Bo nazzola] 
& Gourgoulhon 19961. Furthermore, Jones ( 1975) has pointed out 
that viscous processes can trigger a secular instability, which drives 
the axis of symmetry of the prolate star into the plane perpendicular 
to the angular momentum vector transforming it into a bar-shaped 
rotating source of gravitational radiation (Cutler 2002 , Stell a et aL| 
|2005[ >. From the astrophysical point of view, both the deformation 
of the stellar surface as well as the distortion of the matter dis- 
tribution are relevant and can be measured by appropriate quanti- 
ties, namely the surface deformation (or apparent oblateness) and 
the quadrupole distortion. Previous studies of relativistic models 
of stars with a toroidal magnetic field have provided a broad sur- 
vey of nonrotating and rotating models for varying masses, radii 
and magnetic fluxes. The focus of this work is a complementary 
one to earlier studies: in order to obtain a comprehensive picture 
of the impact of a toroidal magnetic field on relativistic stars, we 
exclusively study models of fixed baryon mass corresponding to 
a gravitational mass of M = 1.400 Mq in the unmagnetised and 
nonrotating case. Although we expect neutrons stars to come in a 
narrow but nonzero range of masses, restricting to a single value 
of the gravitational mass has the important advantage that we can 
explore with unprecedented precision both the deformations intro- 
duced by magnetic fields and those introduced by rapid rotation. 
As we will comment later on, our increased accuracy has allowed 
us also to discover novel results and equilibrium configurations. 

Our neutron stars are modelled assuming the matter to be a 
perfect fluid well described by a single-parameter equation of state 
(EOS), and as having infinite conductivity, as required by the ideal 
magneto-hydrodynamics (MHD) limit. We do not consider multi- 
fluid models, which would allow for the stratification of neutron- 
star matter, nor do we treat the protons in the interior as a super- 
conducting fluid, which would greatly alter the magnetic proper- 
ties of associated equilibrium models. Nevertheless, important re- 
sults obtained for extremely magnetised models with field strengths 
of the order of 10 17 G can be readily extended to configurations 
of much lower and more realistic field strengths of the order of 
10 13 G. More specifically, in addition to a standard y = 2 polytropic 
EOS, we have considered a sample of seven realistic EOSs result- 
ing from calculations of cold catalysed dense matter, namely the 
APR EOS (|Akmal et al.[T998l>, the BBB2 EOS JBaldo et al.| 1997), 



the BN1H1 EOS (Balberg & Gal|1997J, the BPAL12 EOS jBorrT 


|baci|1995|>, the FPS EOS ( Pandharipande & Ravenhall|1989 1, the 


GNH3EOS < 


Glendenning|1985) and the SLy4 EOS (Douchin & 


|Haensel||2001 


i. Altogether, this set of EOSs spans a wide range 



of physical properties and should cover any realistic description 
of neutron stars. For the Pol2 EOS, we have explored systemati- 
cally the space of solutions and computed the corresponding sur- 
face deformation and quadrupole distortion. In addition, we have 
also computed the distortion coefficients, which allow one to com- 
pute the deformation of neutron stars up to large magnetizations 
and rotation rates through a simple algebraic expression, following 
the procedure devised in Cutler (2002). 

The plan of the article is the following: in Sect. [2] we give an 
overview of the novel results obtained in this study, in Sect. [3] we 
discuss the theoretical framework on which our approach, whose 
numerical implementation and testing is discussed in Sect. [4] is 
based. Results for static magnetised models are presented in Sect. [5] 
and for rotating magnetised models in Sect. [6] In Sect. [7] we deal 
with the case of moderate magnetic field and rotation and derive 
empirical distortion coefficients before presenting our conclusions 
in Sect. [8] 



2 OVERVIEW OF THE MAIN RESULTS 

Given the complexity of the numerous results found and the risk 
that the most important ones may be lost in the details, in the fol- 
lowing we briefly summarise what we believe are the most salient 
results obtained. We recall that we measure with e s the deformation 
of the surface shape, while we measure with e the (quadrupolar) 
deviation of the matter distribution from a spherically symmetric 
one. Overall, equilibrium models of relativistic stars with a toroidal 
magnetic field exhibit the following properties: 

• nonrotating and magnetised models exhibit a prolate surface 
deformation and a prolate quadrupolar deformation, i.e. e s < 0, 
6 < 0, both of which decrease as the magnetisation parameter A is 
increased. 

• rotating and unmagnetised models exhibit an oblate surface 
deformation and an oblate quadrupolar deformation, i.e. e s > 0, 
6 > 0, both of which increase as the rotation frequency is increased. 

• between these limiting cases, neutral lines e s = and 6 = 
divide models into having prolate/oblate surface deformations and 
prolate/oblate internal deformations, respectively. 

• for nonrotating models no upper limit was found to the mag- 
netisation parameter Aq, with stellar models that become increas- 
ingly prolate, but also increasingly extended as the magnetisation 
is increased [see Figs.|6]and|7j. 

• the magnetic pressure associated with the toroidal magnetic 
field also causes an expansion in the outer layers of the star, in par- 
ticular a growth of its equatorial radius. This effect is present also 
for nonrotating models, for which, however, the polar radius grows 
more rapidly than the equatorial one, yielding a prolate surface de- 
formation, i.e. 6 S < 0. However, as the rotation is increased and the 
magnetisation decreased, models can be found which have a prolate 
interior deformation, i.e. e < 0, and an oblate surface deformation, 
i.e. 6 S > 0. 

• for any given angular velocity, the model of maximum mag- 
netisation coincides with the mass-shedding model, i.e. the model 
for which centrifugal and gravitational forces are equal at the equa- 
torial radius. This point is characterized by the appearance of a cusp 
at the equator. 

• for rotating models the magnetic pressure at the equator adds 
to the centrifugal force, favouring the loss of mass. As a result of the 
increase in the equatorial size, the mass-shedding frequency is sys- 
tematically smaller than the corresponding one for unmagnetised 
models. 

• in the space of parameters considered, the average toroidal 
magnetic-field strength (B 2 ) 1 ' 2 is not a monotonic function of the 
intrinsic magnetisation parameter A , and after reaching a max- 
imum value, it gradually decreases. This implies the existence 
of double solutions in certain parts of the space of parameters 
(Q 2 ,<B 2 ». 

• although the space of parameters (SI 2 , (B 2 )) is potentially de- 
generate, the corresponding space of parameters (Or, Al) is not. As 
a result, any stellar model considered is characterised uniquely by 
the values of the angular velocity fl and of the magnetisation pa- 
rameter A . 

Given these results, it is natural to divide our models into three 
classes: (1) models labelled PP for prolate-prolate, for which both 
apparent shape and matter distribution are prolate, i.e. e s < 0, e < 0; 
(2) models labelled PO for prolate-oblate, whose shape is oblate 
whereas their matter distribution is prolate, i.e. e s > 0, e < 0; (3) 
models labelled 00 for oblate-oblate, which appear oblate and also 
exhibit an oblate matter distribution; i.e. e s > 0, e > 0. The latter 
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Figure 1. Schematic diagram showing the space of solutions of equilibrium 
models in the (fl 2 ,(B 2 )) plane. According to the relative strength of the 
magnetic field and of the rotation rate, different combinations of the surface 
deformation, e s , and of the quadrupole deformation, e, are possible. 



class of models had not been found by Kiuchi & Yoshida (2008). 
A schematic picture illustrating the three different classes is shown 
in Fig. [I] Finally, while different EOSs with their different stiffness 
introduce quantitative differences in the behaviour described above, 
they all follow the same qualitative behaviour. 



3 MATHEMATICAL AND NUMERICAL SETUP 
3.1 Basic assumptions 

We assume that the spacetime generated by the rotating star is sta- 
tionary and axisymmetric, with Killing vector fields eo and e 3 as- 
sociated with these symmetries. If, in addition, the spacetime is 
asymptotically flat and there exists an axis where e 3 vanishes, then 
eo and e 3 commute ( Carter 1970 ). This enables us to choose coordi- 
nates (x a ) = (t, r, 9, (j>) with vector fields e a = d/dt and e 3 = d/d(f>. If 
furthermore the total stress-energy tensor T satisfies the circularity 
conditions 

T-e = ae +/3e 3 , (1) 
T ■ e 3 = Aeo + \iez , (2) 

convective currents in the meridional planes of constant (/, <f>) are 
absent by construction. In this case, quasi-isotropic (QI) coordi- 
nates can be adopted and the line element reads 

ds 1 = g ap dx^dx 8 = -N 2 dt 2 + OV sin 2 6 2 (d<f> - A"*dr) 2 

+ l F 2 (dr 2 + r 2 d9 2 ) . (3) 

where N, N*, f, and cD are functions of (r, 6). We further intro- 
duce the Eulerian observer Oo whose four-velocity n is the future- 
directed unit-vector normal to hypersurfaces of constant /. From 
Eq. ([3), we infer n a = (-N, 0, 0, 0). 

The compatibility of the electromagnetic fields with the cir- 
cularity condition was established long ago (Carter 1973) for the 
case in which the electromagnetic-field tensor F is derived from 
a potential 1-form A with components (A,, 0, 0, A^), whereas the 



contravariant components of the electric current vector j read f = 
(/, 0, 0, f). This fact has been exploited for computing fully rela- 
tivistic models of stars with a poloidal electromagnetic field ( |Boc-| 
|quet et aLp 995 ; Cardall et al. 2001). However, it has been shown 
that the case of a toroidal magnetic field satisfies the circularity as- 
sumption, too ( |Oron|2002[|Kiuchi & Yoshida|2008) . In the follow- 
ing, we adopt a vector potential which is orthogonal to the Killing 
vectors eo and e 3 , namely A ■ eo = and A ■ e 3 = 0. The covariant 
components of A then become 



A a = (0,A„4,,0). 



(4) 



This particular form of A ensures, by construction, the absence of 
any poloidal electric or magnetic-field component. The only non- 
vanishing component of the antisymmetric Faraday tensor F a p = 



F, u - 



A a fi is then given by 

dA B dA r 
~dr ~ ~dd ' 



(5) 



For the electric field E and the magnetic field B as measured by 
observer Oo, we then obtain 



F afl rP = (0,0,0,0), 



-In 

2 '/ afiyd 

Osinf? 
V 2 



F*rP 



(6) 
(7) 



N* 



dA r 



dAg 
dr 



,0,0,- 



dA r 



dr 



where n n/i y S is the totally antisymmetric tensor associated with the 
metric g. According to Eqs. |6]( and Q, the combination of a van- 
ishing electric field and of a toroidal magnetic field holds for any 
observer whose 4-velocity u is a linear combination of the two 
Killing vectors e and e 3 . The electromagnetic contribution to the 
stress energy tensor is given by 



T, 



0)3 



If 

An \ 



FtyvF" a - 



-F k aF" A g a p 



(8) 



Following the procedure adopted by Bonazzola et af] J1993} , T 
is split up into the total electromagnetic energy density £, the 
Poynting 3-vector J7, and the electromagnetic stress 3-tensor S as 
measured by the Eulerian observer Oo- With the projection tensor 
h = g + n® n, these quantities are obtained as projections of T 
onto and orthogonal to n, namely £ = n ■ T ■ n, = —h-T-n, 
and S = h-T-h, Specialised to the present case of no electric field 
and no poloidal magnetic-field components, the electromagnetic 
contribution to the (3+1) matter variables (for an introduction to 
the (3+1) formalism of general relativity, cf. Smarr & York|1978| 
Gourgoulhon 2012 1, namely the total energy density E = n-T n, the 
momentum density 3-vector J = —h-T ■ n, and the stress 3-tcnsor 
S = h ■ T ■ h reads 



1 



8tt \ O r sin j 

J, = 0, 

*S' j- — & , q — £> , <S^ s — ~& ■ 



(9) 

(10) 
(11) 



In particular, all non-diagonal components of S' j are zero, and <S = 
£. The circularity assumption requires that for the Poynting vector 
J 



J,=0, Je = 0, 

and that the electromagnetic stress tensor S satisfies 



S r4 , = 0, 



-0. 



(12) 



(13) 
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3.2 Einstein Equations 

The Einstein equations for the metric tensor g defined by Eq. |3| 
and a general stress-energy tensor T decomposed into the (3+1) 
quantities E, J, and S become 

(T>2 2 sm 2 q 

Av = An y ¥ 2 (E+S) + — — (dN*) 1 - 8vd(y+B), (14) 

A 3 (A^rsin6») = -16^^-— r sin 6 8^8(38- v), (15) 

<I> 2 rs'md 

A 2 [(iVO- l)rsin6>] = 87rA^F 2 (S r r +S e B ) rsinO, (16) 



A 2 f = 8w>P 2 5^ + 



30>V sin 2 ( 
AN 2 



-{dN*) 2 - (dv) 2 . (17) 



where the following abbreviations have been used 
v = \nN, f = ln(A /> P) , /? = ln<D, 

5^ la i__aj_ 

2 ~ dr 2 + rdr + r 2 88 2 ' 

a 2 2 a i a 2 1 a 

A 3 = + - — + — ^r-r + 



dr 2 rdr r 2 d8 2 r 2 tan 9 39' 

A 3 = A 3 - — -L- , 
r 2 sin 2 

a« az? i a« az? 
0a3bs a7a7 + ^aeae' 



(18) 
(19) 

(20) 

(21) 

(22) 



The resulting system of four nonlinear elliptic equations for the 
metric variables N, A'*, "P, and O can be solved iteratively once suit- 
able boundary conditions of asymptotic flatness have been adopted. 
Additional details can be found in Bonazzola et al.|(l993^ . 

3.3 Maxwell equations 

Since the electromagnetic-field tensor F is derived from a potential 
1-form A, the homogeneous Maxwell equations 



F afi;-y + F/3y;a + Fya fi ~ 



(23) 



are satisfied by construction. The inhomogeneous Maxwell equa- 
tions F a Pj3 = An j" allow to express the electric current 4-vector j 
in terms of the Faraday tensor F, where the alternative expression 



4^/=-L(V=gF^) 



(24) 



with -yFg = v f 2 OA' r 4 sin 9 is used. The electromagnetic field tensor 
F has only one non-vanishing contra- and covariant component, 
which can be now expressed in terms of the azimuthal component 
B$ of the magnetic field B, 



vp2 

a Bit, 



(25) 



(26) 



>F 2 (l>r 2 sin<9 ' 

Combining Eqs. < |24| > and |25j, the poloidal components of the elec 
trie current 4-vector j can be written as 

1 d(NBt) 



Anf 
Anf = - 



V 2 <bNr 2 sin 8 89 ' 
1 WSJ 



(27) 
(28) 



tfPQNr 2 sinO dr ' 
The remaining components / and f are zero, as expected. Note 
that the circularity condition, Eqs. <JTJ, |2](, forbids any meridional 
convective current contributing to /, which is ensured by assum- 
ing the charge carriers to be massless. By taking the divergence of 
Eq. d24b, it follows that /* „ = thanks to the antisymmetry of 



F, and this continuity equation for the electric current leads to a 
dependence between / and f, namely 



V-7 



9 , I .r, 3 , I -H, 



0. 



(29) 



which however is trivially fulfilled by virtue of Eqs. ( |27} and l |28| > 
which express / and f as functions of a single quantity, namely, 
NB,p, without any restriction from Eq. {29} . A useful consequence 
of Eqs. < |27| > and \2S\ is that the flow lines of the electric current 
coincide with the isocontours of NB$, which allows one to visualise 
the current distribution without actually computing j. 

3.4 Equation of motion and condition on the Lorentz force 

In the case of a perfect fluid the stress energy tensor T takes the 
form 



T = (e + p)u u + p g , 



(30) 



where e is the energy density and p the pressure as measured by 
the fluid comoving observer 0\ with 4-velocity u. Since we assume 
the absence of meridional currents, u can be written as a linear 
combination of the Killing vectors «o and e 3 , namely 



■ u eo + ire 3 . 



(31) 



Introducing the Lorentz factor r = -n ■ u = Nit linking observers 
Oo and Oi and the fluid coordinate angular velocity fi = u^/it, the 
physical fluid velocity U and the Lorentz factor T can be expressed 
respectively as 

° rSm6 V^), (32) 



U 



N 



r = (i - u 2 ) m . 



(33) 



Taking into account the contribution of the magnetic field to the 
(3+1) matter variables according to Eqs. |9"|h|1 l|l, the following 
expressions are obtained for a perfect fluid endowed with a toroidal 
magnetic field, 

,2 



E = T 2 (e + p) - p + 



1 



&n \ Or sin ( 
y = r 2 (e + p)cDrsin0t/, 



S'r 



P + 



B. 



p + T 2 (e + p)U l - 



1 

87r \ <br sing 

1 / B. 



4> 



Sn \ <$r sin 9 



(34) 
(35) 
(36) 

(37) 



All other components are zero, and 5 = 3p + T 2 (e + p)U 2 + 
l/(8^)(B^/(Orsine)) 2 . The projection of V • T = onto spatial 
hypersurfaces of constant t provides the equation of motion 



1 



1 



Fmf 



dp dv a „ „ 

+ r :(lnT) rm 

e + p ox' ox' ox' e + p 

where F is defined as 

F = u^u = r — Ur sinO . 



,an 

dx 1 



N 



(38) 



(39) 



Assuming a single-parameter EOS, e = e(n) and p = p(n), with n 
the baryon number density, we introduce the fluid log-enthalpy 



H = In 



e + p 



(40) 



with a mean baryon rest mass of nib = 1.66 x 10~ 27 kg. Using 
Eqs. l |25[H|281 >, l |40[ >, and discarding the case of differential rota- 
tion, Eq. l|38[> becomes 
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0.(41) 



d / 1 W 1 \ d(NB,) 2 

— {H + v-\nT) + \ =- - V- 

dx' \e + P)\ 8/rO 2 r 2 sin 2 ) dx< 

Equation l |41[ > can only be satisfied if the Lorentz force term is de- 
rived from a scalar function M(r, 0) as well, thus we require 



1 



1 



\ d(NB 4> ) 2 



\e + pl\8n$> 2 N 2 r 2 sin 2 



dM 
57 



(42) 



Using the Schwartz theorem, the integrability condition of Eq. {A2\ 
can be written as 



dG d(NBt) dG d(NB?) 



80 dr 



dr 80 



0, 



(43) 



where G = (e + p)(b 2 N 2 r 2 sin' 0. In other words, Eq. 1 43 i states 
that the Jacobian of the coordinate transform (r, 0) — > (G, NB$) is 
zero, and hence there exists a scalar function © : R — > R such 
that &(G,NB^,) = 0. Two different cases can be distinguished: (;') If 
d®ld(NB$) = 0, then does not depend on NB$, and the relation 
&(G, NBj>) = implies that G is constant, which is equivalent to 
the absence of any matter and G = 0. (ii) If 8®ld(NB^) + 0, then 
the implicit-function theorem enables us to write NB# = NB^(G). 
Discarding the case without matter, we retain possibility (ii) and 
conclude that 



NB,f, = g((e + pWN-r sin 2 



(44) 



with g being an arbitrary scalar function. However, the regularity 
properties of an axisymmetric vector field in the case of spatial 
spherical coordinates (r, 0, (j>) impose some further restrictions on 
g. For the covariant comp onent U$ of an azimu thal vector field 
U = U*e 3 , it can be shown ( |Bardeen & Piran|1983"} that 



U^(r,0) = r 2 sin 2 0m(r,0), 



(45) 



where m is an arbitrary axisymmetric scalar function. Application 
of Eq. ( |45| to NBq allows to conclude that g can be written as 
g(x) = xf(x), with an arbitrary regular scalar function /. The re- 
sulting expression for NB$ after application of Eq. \45\ to Eq. l |44| 
reads 

NBp = (e+p) <b 2 N 2 r 2 sin 2 6f((e+p) <5> 2 N 2 r 2 sin 2 0) . (46) 

By inserting Eq. \Ad\ into Eq. \42) , the gradient of the magnetic 
potential M reads 



™ = LJL [(e+p) & N 2 r 2 sin 2 gf]=R . 

ox' An ox' 



(47) 



In general, M cannot be determined in closed form. However, in- 
spection of Eq. \<\1) reveals that for the sub-case of a monomial 
function f(x) = A m x m (not to be meant as a summation over repeated 
indices), a solution can be written down immediately. In this case, 
Eq. d46) simplifies to 



NB, 



X m ((e+p) $ 2 yV 2 r 2 sin 2 0) m+1 . 



(48) 



and the solution to Eq. l |47[ > is given by an algebraic expression, 
namely 



M = —( £±1 1 ((e+p) 9 2 N 2 r 2 sin 2 0) 2 '" +l 
4n \ 2m + 1 1 



(49) 



The simplest function / is obtained for in = and corresponds to a 
constant function of value A . In this case, Eq. (46l reduces to 



NB,p = Ao (e+p) <$> 2 N 2 r 2 sin 2 , 



and thus 



fi = A (e+p) <$> 2 Nr 2 sin 2 



(50) 



(51) 



supplemented by Eq. {49}, which yields the simplified expression 

Al 



M= -r(e+p) $> 2 N 2 r 2 sin 7 
An 



(52) 



For all models computed in this study, the magnetic-field com- 
ponent fi^ is given by Eq. 1 5 1 1 and the magnetic potential M by 



Eq. 1 52 1. Because M vanishes on the axis of symmetry, M c = 0, 
and the integral of Eq. l |41[ > reads 



H + v-lnT + M = H c + v c 



(53) 



Kiuchi & 



where H c and v c are the central values of H and v respectively. Note 
that a general prescription for the determination of the magnetic 
field 6^ and of the magnetic potential M can be found in 
|Yoshida|p008) and |Gourgoulhon et al.|p0TT) . 



3.5 Perfect-conductor relation 

According to Ohm's law and assuming an infinite conductivity for 
the neutron-star matter, the electric field E' as measured by the fluid 
comoving observer O \ has to vanish, namely 



(54) 



E' a = F afs tf = (0, 0,0,0). 

In the case of a purely poloidal magnetic field ( Bocquet et al.|1995] >, 
the perfect-conductor relation Eq. |54j is nontrivial as it establishes 
a dependence between A, and and forces a dependence of the 
angular velocity on A#, i.e. CI = fl(A<j). In the case of a purely 
toroidal magnetic field, on the other hand, the perfect-conductor 
condition is trivially satisfied and no condition is set on the rota- 
tion law, so that differentially-rotating models can be built as in the 
unmagnetised case. 

It is thus possible to consider differentially rotating configura- 
tions following the procedure for the unmagnetised case. The mag- 
netic field B' as measured by the fluid comoving observer Oi now 
becomes 



B' 



l. 



2 ^laPyfi 

T<5> sin0 

VJ/2 

Note that B' + Q.B' 



(55) 



dAo 
dr 



8Ar 

~80 



,0,0,- 



8A r 
~80 



dAo 
dr 



and that B' = TB A 



3.6 Equation of state 

To solve the system of equations introduced in the previous Section, 
we need a prescription relating the energy density e and pressure p 
to the baryon number density n. The simplest of such relations is 
offered by the poly tropic EOS, in which case the following identi- 
ties hold 



e(n) = m^n + k 
p(n) = Kn y , 



7-1 



(56) 
(57) 



where k is the polytropic constant and y the adiabatic index. Com- 
bining Eq. \4Q) and Eqs. l |56[ >, l |57[ >, the log-enthalpy H and n can 
be expressed as functions of each other, namely, as 



(58) 



H(n) = In 1 + — —?—n y - x 
\ m h y-\ 

n(H) = (izli?-^.!)). (59) 

Hereafter, we will assume y = 2 and refer to this EOS as to 
Pol2. In addition to the Pol2 EOS, we have considered a sample of 
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seven single-constituent one-parameter EOS treating the neutron- 
star matter as a perfect fluid and derived from very different models 
of the ground state of cold (zero temperature) dense matter in this 
work as they are provided by the LORENE library. 

More specifically, we have considered the APR EOS ( |Akmal| 
|et al.|1998) , the BBB2 EOS ( |Baldo et al.|1997> , the BN1H1 EOS 
l |Balberg & Gal|1997|>, the BPAL12 EOS pombaci|1993) >, the FPS 
EOS (Pand haripande & Ravenhall|1989| >, the Sly4 EO S (|Douchin| 
|& HaenseipOOl) , and the GNH3 EOS ( |Glendenning]|1985^ . All 
of these realistic EOSs are provided in tabulated form, which re- 
quires the interpolation of listed thermodynamical quantities n, e, 
and p between contiguous sampling points. Thermodynamical con- 
sistency of the interpolated values is ensured by an interpolation 
procedure based on Hermite polynomials introduced by |Swesty| 
( 19961, which is crucial for the accuracy of the resulting numeri- 
cal models < |Salgado et al.|1994J . 

The properties of the spherical nonrotating and unmagnetised 
reference models with a gravitational mass of M = 1.400 Mq can 
be found in Table[5] The softest EOS of this sample is the BPAL12 
EOS, which yields the smallest circumferential radius of ^ cilc = 
10.06 km, whereas the stiffest one is the GNH3 EOS, which yields 
the largest circumferential radius of R C1TC = 14.20 km. 



3.7 Global quantities 

A number of global quantities which characterise the numerical 
models presented in this study can be computed and will be needed 
in the course of this investigation. These are: the gravitational mass 



where r c , r p are the equatorial and polar coordinate radii, respec- 



M 



E + S + — Ntjjr sine dr 60 1 



the total angular momentum 
J= f Y^JtSsmOdrdOdt, 
the total magnetic energy 

^QSfsmOdrdOdcb 



the total kinetic energy 
1 

T = -a/, 

2 

and the gravitational binding energy 
W = M - T - M p - J£ , 

where the total internal energy M p is defined as 



M, 



Ter s'mOdrdQt 



while the total baryon mass of the star is given by 



M h 



Tnr sin0drd6>c 



(60) 



(61) 



(62) 



(63) 



(64) 



(65) 



(66) 



Other important quantities that are more closely related to the 
deformation of the star are the circumferential radius R Q \ m , which is 
the equatorial radius as measured by the observer Og and is related 
to the equatorial coordinate radius R according to 



R Ac =<I>(R,n/2)R, 

the surface deformation (or apparent oblateness) 
e s = r c /r p - 1 , 



(67) 
(68) 



tively, and the quadrupole distortion e of the star ( jBonazzola &| 
|Gourgoulhon|1996| > 



3 

2~r 



(69) 



where J'- is the quadrupole moment measured in some asymp- 
totically Cartesian mass-centred coordinate system jThorne|1980| >, 
while the moment of inertia / = /.. is defined as / = J /SI. We stress 
that no moment of inertia other than the latter can be defined in a 
meaningful way for axisymmetric rotating bodies and that the rota- 
tion must be rigid and about the axis of symmetry. Furthermore, we 
note that -(3/2)J? z , differs from the standard quadrupole moment 
Q and that in the case of QI coordinates, the latter can be read off 
from the asymptotic expansion of In N according to 



M 



M 2 
12? 



i-- 1 



(70) 



Thome's quadrupole moment is the relevant quantity when the 
gravitational-wave emission from a distorted star has to be deter- 
mined, and the corresponding values are extracted from the asymp- 
totic expansion of certain components of the metric tensor g in 
QI coordinates following the procedure presented in Bona zzoia &| 
|Gourgoulho"nl(T996) . 



4 NUMERICAL IMPLEMENTATION AND TESTS 
4.1 Numerical scheme 

The analytic scheme presented in Sect. [3] has been implemented 
by extending an existing multi-domain, surface-adaptive spectral 
method for unmagnetised relativistic stars presented in |Bonazzola 
et al.|jl998) ; |Gourgoulhon et al.|fi999} and part of the L0REN 



C++ class library for numerical relativity. We refer to Bonazzoia 



|et al.|{T998| l; |Gourgoulhon et al.| ( |1999| > for a detailed description 
of the numerical method and of the tests carried out. The numerical 
solution is computed by iteration, starting from crude initial con- 
ditions of a parabolic log-enthalpy profile, no rotation, flat space 
and no magnetic field. Convergence is monitored by computing 
\\H„ - #„-ill/||#„_i|| where ||if„_i|| is the sum of the absolute values 
of the log-enthalpy H over all collocation points at iteration step 
n — 1. A convergent solution is assumed to be found when the nor- 
malized log-enthalpy residual decreases below a prescribed thresh- 
old, which is usually chosen to be of the order of 10~ 8 or smaller. 
The computational domain is divided into a spherical nucleus con- 
taining the star and two shells covering the exterior, where the last 
one maps onto a finite interval the whole exterior space, from a 
certain radius R up to spatial infinity. The default number of col- 
location points used in our study has been n e = 17 in the angular 
direction and n R = (33, 33, 17) in the radial direction, where the dif- 
ferent numbers refer to the nucleus and the two shells, respectively. 
However, depending on the circumstances, the numerical resolu- 
tion has been increased when necessary, i.e. up to n e = 129 and 
n R = (257, 257, 129) in the case of huge models with a matter dis- 
tribution strongly deviating from a spherical one [cf. Tableland 
Fig.§. 
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2 


5/3 


3/2 


4/3 


9/7 


Sood and Trehan (1972) 
This work 


(e/h) 


-0.21377 
-0.21376 


-0.14284 
-0.14218 


-0.09292 
-0.09292 


-0.03374 
-0.03374 


-0.01724 
-0.01724 



Table 1. Comparison of results for the normalized oblateness e s //i computed in the Newtonian limit for various polytropic EOS with results from |Sood &| 
Trehan 19721 (values for e/h from table I where e = (r e — r p )/ro and ro is the radius of the unperturbed model). The perturbation parameter h = A^/An 1 is a 



measure of the strength of the magnetic field. 



4.2 Comparison with Newtonian results 

The results of the numerical scheme in the Newtonian limit have 
been compared with those of an earlier investigation in the New- 
tonian regime. [Sood & Treh an ( 19721 have in fact computed lin- 
ear perturbations of polytropic stars induced by a toroidal magnetic 
field corresponding to the Newtonian limit of Eq. \5\) , namely 

— t —=A a prsm8, (71) 
rsinfc* 

where p denotes the Newtonian mass density. In Table[T| we show 
normalized values of the resulting oblateness for different adiabatic 
exponents produced with the present code and the corresponding 
results from Sood & Trehan ( 1972). We note that |Sood & T rehan 
( |1972) > measured the deformation as e s = (r c - r p )/ro where ro is 
the radius of the unperturbed model. Values computed according 
to this definition coincide with e s = r e /r p - 1 within the rounding 
error of tabulated values. Clearly, for the perturbation parameter 
h = A^/4n 2 , the respective values of e s /h agree to within 10" 5 , 
exceptfory = 5/3, for which the difference is about 10 -3 . However, 
for y = 5/3, certain functions intervening in the solution of the 
perturbed Lane-Emden equation become singular at the surface of 
the star ( Das & Tandon 1977 ), suggesting a larger numerical error 
for the y = 5/3 case as computed by[Sood & Trehan ( 1972| >. 

4.3 Virial identities 

In order to monitor the global error of the numerical models com- 
puted with the new code we have calculated the two general- 
relativistic virial identities GRV2 (Bonazzola 1973, Bonazzola & 
|Gourgoulhon|1994) and GRV3 ( |Gourgoulhon & Bonazzola|1994| >. 
We recall that these identities have to be fulfilled by any solution 
to the Einstein equations Eqs. j!4|l-jl7|l and are not enforced dur- 
ing the iterative procedure that leads to the solution. The GRV2 
identity, in particular, has been shown to be directly related to the 
global error of a numerical solution I Bonaz zola et al.| 1 993 ), while 
the GRV3 identity is an extension of the virial theorem of New- 
tonian physics to general relativity [see Nozaw a~et al.| ( [l998) for 
details of the practical implementation]. For strongly-magnetised 
nonrotating models based on a polytropic EOS with y = 2 and the 
moderate number of collocation points specified in Sect. |4. 1| the 
violations to the identities GRV2 and GRV3 have been found to 
be the order of 10 6 ore better. For the huge models listed in Ta- 
ble]^ corresponding values are of the order of 1CT 4 which required 
the number of collocation points to be multiplied by a factor eight 
compared to their default values. In the case of rapid rotation, on 
the other hand, the numerical grid is not perfectly adapted to the 
stellar surface, and the spectral approximation suffers from Gibbs 
phenomena, as reported in Bonazzola et al. ( 1993 ). However, the 
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corresponding values of the identities GRV2 and GRV3 are still of 
the order of 10~ 5 . For realistic EOS, values of GRV2 and GRV3 are 
of the order of 1CT 4 or better. 



5 NONROTATING MAGNETISED MODELS 

Nonrotating models of neutron stars with a toroidal magnetic field 
always generate static spacetimes for which the time Killing vector 
«o is hypersurface-orthogonal. Note that this assumption does not 
always hold in the case of a poloidal magnetic field because the 
additional presence of an electric field gives rise to an azimuthal 
Poynting vector [cf. Bona zzola et al.|{l993[ >], which introduces a 
nonzero angular momentum even when the star does not rotate. As 
a consequence, the shift vector component does not vanish. 

Although we are interested in the construction of rapidly ro- 
tating models (that we will present in Sect.[6](, nonrotating configu- 
rations allow us to investigate the effects of the magnetic field with- 
out the influence of rotationally-induced effects. As a representative 
example, Fig. [2] shows isocontours of the magnetic-field strength 
and baryon number density for the nonrotating model PP-Pol2 
built with the Pol2 EOS and for the maximum field strength (the 
physical properties of model PP-Pol2 are listed in Table|4j. The 
unmagnetised reference model with a circumferential radius of 
R C1IC = 12.00km is symbolised by a gray disc. In agreement with 
the simple structure function defined in Eq. ( |51fr , the toroidal mag- 
netic field shown in the left panel vanishes on the axis of symmetry, 
reaches a maximum value of S m;lx = 7.408x 10 17 G in the equatorial 
plane deep inside the star, and decreases towards the surface of the 
star, where it vanishes, so that the magnetic field is fully contained 
inside the star. The forces exerted by the magnetic field can be vi- 
sualised through the magnetic potential M, which shows a distribu- 
tion similar to that of the magnetic-field strength and that, for this 
reason, we do not report here. It should be noted that the magnetic 
potential M is repulsive, since it is positive everywhere by construc- 
tion. This is to be contrasted with the magnetic potential associated 
with the purely poloidal magnetic field adopted in Bocqu eTet al.| 
( [19951 , th at 

was instead attractive. As a result, the corresponding 
forces are the opposite, despite the isocontours appear to be very 
similar. Note also that in the present case of a purely toroidal mag- 
netic field defined according to Eq. |5l| , the isocontours of M co- 
incide with the flow lines of the electric current j, which are nested 
loops in the meridional plane. Finally, we note that M vanishes at 
the surface of the star, implying that the latter coincides with an 
isosurface of the gravitational potential v, as required by Eq. l |53[ ) 
whenlnT = and M = 0. 

The forces that can be derived from M vanish on the axis of 
symmetry and reach their maximum in the equatorial plane. Be- 
tween the centre of the star and the maximum of M, they are di- 
rected inwards, inducing an approximately cylindrical compression 
of the central region of the star, which responds by a prolate defor- 
mation visible in the right panel of Fig. [2] In the outer layers of 
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Figure 2. Isocontours of magnetic-field strength (left panel) and baryon number density (right panel) in the (x,z) plane of model PP-Pol2 of a nonrotating 
star with a gravitational mass of M = 1.400 Mq and a circumferential radius of R c j rc = 12.00 km in the unmagnetised case built with a polytropic EOS with 
y = 2 for a maximum average magnetic-field strength of (B 2 ) 1 ' 2 = 2.917 X 10 17 G. The gray disc indicates the dimensions of the unmagnetised reference 
model. Physical properties of model PP-Pol2 are listed in Table|4] 





Figure 3. Comparison of the gravitational mass M, circumferential radius R c j rc , central baryon number density n bf , and mean deformation e along with the 
maximum magnetic-field strength B max attained inside the star as a function of the magnetisation parameter Aq for a Pol2 EOS model with a baryon mass of 
M b = 1 .680 Mq and a circumferential radius of K c j rc = 14.30 km in the unmagnetised case with results of Kiuchi & Yoshida 1 2008 1. The relative variation of 
some quantity^ is defined as A^(B max , M b ) = Cr(B max , M b ) ~x(0, M h J)/x (0, M b ). 
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Figure 4. Isocontours of magnetic-field strength (left panel) and baryon number density (right panel) in the (x, z) plane of model PP-APR of a nonrotating 
star with a gravitational mass of M = IAOOMq and a circumferential radius of i? c j rc = 11.34km in the unmagnetised case built with the APR EOS for a 
maximum average magnetic-field strength of (B 2 ) 1 ' 2 = 3.597 X 10 17 G. The gray disc indicates the dimensions of the unmagnetised reference model. The 
physical properties of model PP-APR are listed in Table|4] 



the star, on the other hand, the magnetic forces are directed radi- 
ally outward, pushing them away from the centre, which results in 
a growth of the dimensions of the star whose circumferential radius 
has increased now to a value of R ciTC = 14.34 km. 

Both the apparent shape and the matter distribution are prolate, 
corresponding to negative values of the surface deformation and of 
the quadrupole distortion, namely e s = -0.0806 and e = -0.1986, 
respectively. The value of e s is significantly smaller than that of e, 
and this due to the fact that terms of higher order in the multipole 
expansion of the gravitational potential v fall off rapidly while, at 
the same time, the surface of the star must coincide with an isosur- 
face of v. This interpretation is supported by the fact that this dif- 
ference is significantly smaller for small perturbations around this 
model, as confirmed by comparing values of respective distortion 
coefficients b B and c B reported in Sect. [7] 

The strength of the magnetic field is controlled by the mag- 
netisation parameter A [cf. Eq. \52\\, and the physical quantities 
associated with our models can be parametrized accordingly, in 
particular the maximum magnetic-field strength B max . In Fig. [3[ 
we compare our results with those reported in fig. 6 of Kiuchi & 
|YoshidaH 2008F| considering in particular the variation of the grav- 
itational mass M, of the circumferential radius R CIIC , of the central 
baryon number density n b c , and of the mean deformation e. The 
comparison is made as a function of the maximum magnetic-field 
strength B max and for a nonrotating reference model built with the 
Pol2 EOS and having a baryon mass of M b = 1.680 Mq and a 
circumferential radius of i? c i rc = 14.30 km in the unmagnetised 
case. For each of these quantities we define the difference as 
A y(£ max ,M b ) = \x(B max ,M b ) - x(0,M b )]/x(0,M b ). We note that 
in |Kiuchi & Yos hida (2008) the mean deformation is defined as 
e = (l zz — 1 XX )/I ZZ , where l xx ,l zz are the momenta of inertia rela- 
tive to the corresponding axes [cf. Eq. (3.12) of Kiuchi & Yoshida 



(2008 )]; this measure is similar but different from our measure of 
the quadrupole distortion ^] 

The most prominent feature of this comparison is that B max 
is not a monotonic function of Aq, thus being responsible for the 
presence of a turning point located at a magnetic-field strength of 
Smax = 6.141 x 10 17 G. While our results (blue solid lines) agree 
qualitatively with those of Kiuchi & Yoshida ( 2008 1 (red dashed 
lines), we find significant quantitative differences for all of the 
quantities considered, the largest one being an increase in n bjC of 
7% against the smaller 3% increase found by Kiuchi & Yoshida 
(2008). Another important difference is in the values of the maxi- 
mum magnetic-field strength, which is B max = 6.141 x 10 17 G for 
us and B max = 5.503 x 10 17 G for |Kiuchi & Yoshida| ( |2008| >. The 
differences between the two calculations are smaller when compar- 
ing the mean deformations e, but only for moderate values of B max . 
Furthermore, the measure of e is considerably different from the 
quadrupole deformation e (green dotted line), suggesting that the 
definition of e is not suitable for estimating the gravitational-wave 
emission of a distorted star because it overestimates corresponding 
values by about a factor of two for this reference model. Similar 
differences in the values of e and e have been found also for the 
unmagnetised model with a baryon mass of M b = 1.780 Mq ro- 
tating at fi = 3.230 x 10 3 s _I from table IV of |Kiuchi & Yoshida 
(2008); in this case, the respective values of the mean deforma- 
tion of e = 0.07764 and e = 0.07462, differ by about 4%, the 
quadrupole distortion for the same model is e = 0.03870 which is 
again only half of the value of e. 

In addition to the polytropic Pol2 EOS, we have computed 
nonrotating models for the sample of realistic EOSs discussed in 
Sect. |3.6| and for reference models having a gravitational mass of 
M = 1 .400 Mq in the unmagnetised case (basic physical properties 
of the unmagnetised reference models are collected in Table[5]l. Be- 



2 The data from fig. 6 of Kiuchi & Yoshida (2008j has been read off with 
the utility g3data, see https://github.com/pn22S8/g3data. 



3 More specifically, the definition of & is valid only in a Newtonian frame- 
work, where (/. ; - I XX )II~ = (J - Ija)/^ = 2Newt/A and QNewt is the New- 
tonian quadrupole moment. 
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cause we cannot show equilibrium models for all of these models, 
we have decided to use as alternative reference EOS the APR one 
and to complement model PP-Pol2 at the maximum field-strength 
limit, with the corresponding model PP-APR, which nicely illus- 
trates the impact of a prototypical realistic EOS on the resulting 
equilibrium model. 

Figure [4] shows isocontours of magnetic-field strength and 
baryon number density for model PP-APR at the maximum field- 
strength limit, whose physical properties are collected in Table [4] 
(again, the dimensions of the unmagnetised reference model with 
a circumferential radius of R circ = 11.34 km are indicated by a 
gray disc). The maximum magnetic-field strength attains a value 
of Bmx = 8.046 x 10 17 G in the equatorial plane. Because the 
APR EOS is stiffer than the Pol2 EOS, the matter distribution 
appears less condensed and the peaks of the magnetic field have 
moved slightly outwards. Moreover, both the surface deformation 
6, = -0.1 176 and the quadrupole distortion e s = -0.3045 are larger 
than the corresponding values of model PP-Pol2 by ~ 50%. In 
contrast with model PP-Pol2, a substantial fraction of the stel- 
lar interior below the surface appears to be field-free. This is easy 
to understand: according to Eq. l |51[ >, in fact, the amplitude of the 
toroidal magnetic field is proportional to (e + p), so that the pres- 
ence of a low-density crust as in the APR EOS suppresses the pres- 
ence of a toroidal magnetic field in the outer layers of the star. The 
contour of the unmagnetised reference model not only emphasises 
the significant prolate deformation of the maximum field strength 
model, but also the important increase in the dimensions of the star. 

Overall, the nonrotating models built with realistic EOSs show 
a behaviour which is qualitatively similar to that of the polytropic 
Pol2 EOS as shown in Table [5] and, for increasing magnetisation, 
all realistic EOS exhibit a maximum value of the magnetic-field 
strength (B 2 ) 1 /2 , beyond which it then decreases. The smallest value 
is obtained for the GNH3 EOS, with <B> = <B 2 > 1/2 = 2.161 x 10 17 
G; the largest one is instead obtained for the BPAL12 EOS, with 
{B 2 ) lj2 = 4.790 x 10 17 G. Relevant data for all maximum field 
strength models is collected in Tableland reveals that the values of 
the maximum magnetic-field strength decrease with increasing cir- 
cumferential radii B c j rc , and which itself is a telltale of the stiffness 
of the EOS. As a result, the soft BPAL12 EOS with a circumferen- 
tial radius of only B cin . = 1 1.83 km yields the model with the high- 
est maximum magnetic-field strength, while the stiff GNH3 EOS 
the model with the lowest one. 

We note that in Fig. [3] we have employed the peak magnetic- 
field strength B max and not the average magnetic-field strength 
(B 2 ) 1 ' 2 , which reaches its turning point already at a lower magneti- 
sation level. For this reason, all models based on a realistic EOS are 
located in a region where An b >c > 0, as can be seen by a compari- 
son with the related values of the unmagnetised models compiled in 
Table [5] The maximum values of the peak magnetic-field strength 
B max are sensibly larger than those of the average magnetic-field 
strength (B 2 ) 1 ' 2 , attaining their maximum value for the model built 
with the FPS EOS with B max = 1.210 x 10 18 G. 

The slight increase of the gravitational mass is the result of dif- 
ferent contributions, which can be illustrated by the example of the 
PP-Pol2 model. More specifically, the increase by 0.02682 Mq in 
the gravitational mass of the maximum field-strength model with 
respect to the unmagnetised model. This increase is the sum of a 
positive contribution of 0.03036 Mq due to the total magnetic en- 
ergy, of a negative contribution of -0.02160 Mq due to the internal 
energy lost because of the increased volume, and of a positive con- 
tribution of 0.01806 Mq by which the magnetised model is less 
gravitationally bound than its unmagnetised counterpart. Further- 



more, all maximum field-strength models exhibit smaller momenta 
of inertia than in the unmagnetised case because of the lateral com- 
pression of the stellar core. Finally, we note that after reaching a 
minimum value, the momenta of inertia increase continuously with 
the magnetisation, as the growth of the dimensions of the star be- 
comes significant. 

The values of the surface deformation e s and of the quadrupole 
distortion e obtained for our sample of maximum field-strength 
models also show a noticeable dependence on the stiffness of the 
EOS. In particular, the value of e decreases (its absolute value in- 
creasing) with the increase of the circumferential radius B circ . The 
Pol2 EOS model appears somehow at variance with this trend, but 
we believe this is caused by our choice of a comparatively large cir- 
cumferential radius of B cilc = 12.00 km in the unmagnetised case, 
and which does not reflect the rather soft character of this EOS. 
Lower and upper bounds for the quadrupole distortion are set by the 
two extremal EOS of our sample: the model built with the soft EOS 
BPAL12 with a radius of less than 12.00 km exhibits the smallest 
quadrupole distortion despite of the largest ratio of magnetic energy 
to binding energy, whereas the large GNH3 model with a radius of 
17.06 km shows the largest quadrupole distortion for the lowest ra- 
tio of magnetic energy to binding energy of the whole sample. 

Much of what discussed so far is summarised in Fig. [5] where 
we present the dependence of the surface deformation e s and of 
the quadrupole distortion e on the average magnetic-field strength 
<B 2 > 1/2 , for the representative Pol2 EOS and the APR EOSs. Note 
that up to (B 2 ) 112 ^ 2.0 x 10 17 G, the behaviour of both distor- 
tions is essentially the same for the two EOSs and an almost per- 
fectly linear functiorj^] of (B 2 ). This is not particularly surprising 
and indeed earlier Newtonian studies (Wentzel I960, Ostriker & 
|Gunn|1 969) have suggested to parametrize the quadrupolar distor- 
tion eM ewt induced in a self-gravitating incompressible fluid by a 
toroidal magnetic field as a simple function of the ratio of the mag- 
netic energy to the binding energy, and hence as a function which 
is quadratic in the magnetic-field strength. Figure[5]suggests there- 
fore that this behaviour is preserved also in general relativity and 
up to very large magnetic fields, thus offering the possibility of ex- 
pressing the magnetic-induced deformation in terms of a simple 
algebraic expression with coefficients which will correct the previ- 
ously known Newtonian ones. A more detailed discussion of this 
point is presented in Sect.|7]and in Appendix [A] 

We conclude this Section by commenting on a novel and par- 
ticularly interesting result of our analysis of nonrotating configura- 
tions, namely, that no physical or numerical limit was found for the 
magnetisation level. More specifically, after suitably increasing the 
number of collocation points according to the desired level of ac- 
curacy, we were able to find solutions with very large values of the 
surface and quadrupolar deformations. As an example, we could 
obtain convergent solutions of the nonrotating reference model of 
|Kiuchi & Yoshida|p008l > with a baryon mass of M b = 1.680 Mq 
and a circumferential radius of B circ = 14.30 km in the unmagne- 
tised case, now strongly magnetised and attaining a circumferential 
radius of R c „ c = 101.5 km and a related value of AB c j lc = 6.098. 
While the average magnetic-field strength (B 2 } l/2 = 0.1461 x 10 17 
G is about an order of magnitude smaller than the maximum value 
of (B 2 ) 1 ' 2 = 2.241 x 10 17 G reached at a circumferential radius of 
B circ = 17.09 km, the ratio of magnetic energy to binding energy is 



This is not straightforward to deduce from our Fig. |5j which reports e s 
and e as a function of (B). However, it is very apparent when plotting the 
distortions as a function (B 2 ) (not reported here for compactness). 
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(B 2 )" 2 [10 17 G] <B 2 > 1/2 [10 17 G] 

Figure 5. Surface deformation e s (left panel) and quadrupole distortion £b (right panel) for nonrotating models built with the Pol2 EOS and the APR EOS as a 
function of the root mean square magnetic-field strength (B 2 ) 1 ' 2 from the unmagnetised limit up to the maximum field strength models PP-Pol2 and PP-APR, 
respectively. The physical properties of models PP-Pol2 and PP-APR are listed in Table|4] 



Table 2. Nonrotating neutron star models for large magnetisations of .4? j\W\ » 0.5. The gravitational mass is M = 1 .400 Mq for each EOS respectively in 
the unmagnetised and nonrotating case for which properties are listed in Table[5] B max is the maximum value of the magnetic field, (B 2 ) 1 ' 2 the root mean 
square average of the magnetic-field strength determined over the volume of the star, W\ the absolute value of the ratio of total magnetic energy ^ and 
potential energy W, «b,c the central baryon number density [0.1 fm~ 3 ], M the gravitational mass, B c j lc the circumferential radius, / the moment of inertia, e s 
the surface deformation, e the quadrupole distortion, and GRV2/GRV3 the estimates of the global error of respective models based on the relativistic virial 
identities introduced in Sect. 14. 31 
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much higher with values of JCIW\ = 0.4868 and JK/\W\ = 0.1308 
respectively. The surface deformation of this model reaches a value 
of 6 S = -0.2403 and the quadrupole distortion an impressive value 
of e = -6.946. While clearly unrealistic for ordinary neutron stars, 
young proto-neutron stars may, in principle, attain these magneti- 
zations and possibly these sizes 1 Vil lain et al.|20 04); clearly, this 
represents a possibility that deserves a more careful investigation. 

In order to illustrate the extreme effects of the toroidal mag- 
netic field at high levels of magnetisation, we have computed 
model MM-Pol2 of the Pol2 EOS reference model with a gravi- 
tational mass of M = 1.400 Mq and a circumferential radius of 
R cilc = 12.00 km in the nonrotating and unmagnetised case, now 
inflated to a circumferential radius of R cill . = 100.0 km. Such an 
"extra-large" model is shown in Fig. [6] which, as the previous ones, 
reports the isocontours of magnetic-field strength (left panel) and 
baryon number density (right panel) for an average magnetic-field 
strength of (B 2 )' /2 = 0.1400 x 10 17 G. The physical properties of 
model MM-Pol2 are listed in Table [2] In this case, therefore, the 
surface and quadrupolar deformations reach the extreme values of 
6 S = -0.2609 and of e = -7.803, respectively. 

Although new and somewhat surprising, these results are not 
totally unexpected, and we note that already for a purely poloidal 
magnetic field, Card all et al.| ( |2001[ l doubted the existence of a 
mass-shedding limit in the nonrotating case. This leads to the con- 
clusion that the non-convergence limit encountered by |Kiuchi &| 
Yoshida (20081 at R ciTC = 28.85 km may be due to their numerical 
scheme. Furthermore, to confirm that the extreme distortions re- 
ported for model MM-Pol2 are not a peculiarity of the Pol2 EOS, 
we have found similar equilibria also for model MM-APR, built 



with the APR EOS with a gravitational mass of M = 1.400 Mq 
and a circumferential radius of R cill . = 11.34 km in the unmagne- 
tised case. This is shown in Fig. [7] and exhibits the same spindle- 
shaped matter distribution of the stellar core for a circumferen- 
tial radius of R CKC = 90.26 km and an average magnetic-field 
strength of <B 2 > 1/2 = 0.2313 x 10 17 G. The physical properties 
of model MM-APR are listed in Table [2] but we note here that the 
surface and quadrupolar deformations reach the extreme values of 
6 S = -0.1248 and of e = -4.810, respectively. These values are 
smaller than those of model MM-Pol2 (despite the ratio of mag- 
netic energy to binding energy is similar and ^£ l\W\ ~ 0.5) but 
shows quite clearly that extremely large deformations are a feature 
of nonrotating models, independently of the EOS. 



6 ROTATING MAGNETISED MODELS 

Having investigated in the previous Section nonrotating models, we 
next turn to models that include rotation, which is well known to 
induce deformations of the stars that are the opposite of those dis- 
cussed so far, namely, of introducing an oblateness generically both 
in the surface deformation and in the quadrupole distortion. In par- 
ticular, we are concerned with determining the limits of the space 
of solutions both in terms of the magnetisation and of the rota- 
tion rate. To probe in detail such a space of solutions we consider 
two rotating models corresponding to maximum field-strength con- 
figurations of the Pol2 EOS. The first of these reference models, 
which we refer to as P0-Pol2, has a mean magnetic-field strength 
of (B 2 ) 1 ' 2 = 2.324 x 10 17 G and rotates at an angular velocity of 
£2 = 3.969 x 10 3 s , thus with a moderate ratio of kinetic en- 
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Table 3. The same as in Table[2]but for models computed with different EOSs at the maximum field-strength limit. 
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Figure 6. Isocontours of magnetic-field strength (left panel) and baryon number density (right panel) in the (x, z) plane of model MM-Pol2 of a nonrotating star 
with a gravitational mass of M = 1.400 Mq and a circumferential radius of R c j lc = 12.00 km in the unmagnetised case built with the Pol2 EOS for an average 
magnetic-field strength of (S 2 > 1/2 = 0.1400 X 10 17 G. The physical properties of model MM-Pol2 are listed in Table[2] 




x [km] x [km] 

Figure 7. Isocontours of magnetic-field strength (left panel) and baryon number density (right panel) in the (x, z) plane of model MM-APR of a nonrotating star 
with a gravitational mass of M = 1 .400 Mq and a circumferential radius of R mIC = 1 1 .34 km in the unmagnetised case built with the APR EOS for an average 
magnetic-field strength of (B 2 ) 1 ' 2 = 0.2313 X 10 17 G. The physical properties of model MM-APR are listed in Table^ 
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Figure 8. Isocontours of magnetic-field strength (left panel) and baryon number density (right panel) in the (x,z) plane of the P0-Pol2 and the 00-Pol2 
models of a star built with a polytropic EOS with y = 2 with a gravitational mass of M = 1.400 Mq and a circumferential radius of R t .j rc = 12.00 km in the 
unmagnetised and nonrotating case which is now rotating at £1 = 3.969 X 10 3 s - ' (top) and CI = 5.050 X 10 3 s (bottom) respectively. The gray disc indicates 
the dimensions of the unmagnetised reference model. The physical properties of models P0-Pol2 and 00-Pol2 are listed in Table|4] 



ergy to binding energy of 77|W1 = 0.03200. The second model, 
that we refer to as 00-Pol2, has a mean magnetic-field strength of 
(5 2 ) 1 / 2 = 1.393 x 10 17 G and rotates rapidly at an angular velocity 
of fi = 5.050 x 10 3 s~\ with a ratio T/\W\ = 0.07136 (see Table|4] 
for a complete list of the physical properties). 

The magnetic-field strength and the baryon number density 
of both models are shown in Fig. [8] with the top row referring to 
P0-Pol2 and the bottom one to 00-Pol2. The denomination of 
these two models becomes apparent when considering the corre- 
sponding surface and quadrupole deformations. While in fact the 
qualitative properties of P0-Pol2 seem to resemble those of the 
nonrotating model PP-Pol2, the surface of this model is markedly 
flattened and indeed with a positive oblateness of e s = 0.1521. The 
matter distribution inside the star, however, is still prolately de- 
formed, as confirmed by a negative value of the quadrupole distor- 
tion, namely e = -0.0684. On the other hand, for model 00-Pol2 



both the surface deformation and the quadrupole distortion are pos- 
itive, i.e. e s = 0.5861 and e = 0.0856, even though the isocontours 
of the baryon number density are still prolate toward the centre 
of the star. Interestingly, for model P0-Pol2, for which the mag- 
netic field is still dominating, the ratio of the magnetic energy to 
the binding energy /\W\ exceeds that of the kinetic energy to the 
binding energy T/\ W\, whereas for model 00-Pol2, the opposite is 
true. Hence, while this condition does not hold in all cases, we can 
take the inequality ./£ '/T > 1 as a first approximate criterion for the 
production of a negative quadrupole distortion. A more quantitative 
discussion on this will be presented in Sect. [7] 

In our sampling of the space of parameter we have computed 
a total of more than 900 models of rotating and magnetised equilib- 
rium configurations of the Pol2 EOS, which are uniquely labelled 
by the values of the angular velocity Si and of the magnetisation 
parameter A . The latter extends up to a maximum obtained for a 
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Figure 9. Isocontours of the surface deformation e s for the Pol2 EOS reference model with a gravitational mass of M = 1 .400 Mq and a circumferential radius 
of S c j rc = 12.00 km in the unmagnetised and nonrotating case. The curved line which connects the nonrotating model with a maximum average magnetic-field 
strength of (B 2 ) i/2 = 2.917 X 10 17 G with the unmagnetised model rotating at the mass-shedding limit of CI = 5.205 X 10 3 s 1 separates the lower part of the 
solution space (left panel) where the mean magnetic-field strength increases as a function of the field strength parameter Aq from the upper part of the solution 
space (right panel) where the mean magnetic-field strength decreases as a function of the latter. In the right panel, we have also shown the lower sheet lying 
underneath. 




Figure 10. Isocontours of the quadrupole distortion e for the Pol2 EOS reference model with a gravitational mass of M = 1 .400 Mq and a circumferential 
radius of R c j lc = 12.00 km in the unmagnetised and nonrotating case. The curved line which connects the nonrotating model with a maximum average 
magnetic-field strength of (B 2 ) ,/2 = 2.917 X 10 17 G with the unmagnetised model rotating at the mass-shedding limit of CI = 5.205 X 10 3 s 1 separates the 
lower part of the solution space (left panel) where the mean magnetic-field strength increases as a function of the field strength parameter Aq from the upper 
part of the solution space (right panel) where the mean magnetic-field strength decreases as a function of latter. In the right panel we have also shown the lower 
sheet lying underneath. 



nonrotating model with a radius of R ciTC = 19.45 km and a ratio of 
magnetic energy to binding energy of J?/\W\ = 0.244^] Overall, 
the space of physical solutions is delimited by four boundaries: (1) 
the nonrotating limit with Q. = 0; (2) the unmagnetised limit with 



Note that the "extra-large" models shown in Figs. |6jand|7] are located 
beyond the truncation limit and thus do not belong to what we consider the 
parameter space. 



6 = 0; (3) the (self-imposed) "truncation limit" with respect to the 
magnetisation parameter Aq; (4) the mass-shedding limit, beyond 
which no rotating solution does exist. 

We have already noted that for nonrotating models the mean 
magnetic-field strength (B 2 } 1 ' 2 is not a monotonic function of the 
magnetisation parameter A and that after attaining a maximum 
value of (S 2 ) 1 ' 2 = 2.917 x 10 17 G, it decreases continuously as the 
magnetisation is increased. This behaviour has been found also for 
rotating models, thus implying a non-uniqueness for models in the 
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space (S1 2 ,(B 2 )). Such a degeneracy could be avoided by replac- 
ing the average magnetic-field strength with a quantity that grows 
monotonically with the magnetisation, e.g. the circumferential ra- 
dius Rex,, or the ratio of magnetic energy to binding energy ,/£ l\W\. 
However, because of their fundamental astrophysical importance, 
we have chosen to report the results for the surface deformation e s 
and the quadrupole distortion e in terms of the mean magnetic-field 
strength (B 2 ) 1 ' 2 as the ordering quantity and of the angular veloc- 
ity SI. The presence of this degeneracy implies that when evaluated 
in the space of parameters (CI 2 , (B 2 )), the distortions e s and e will 
select a two-dimensional surface which can be split up along the 
turning points of maximum magnetic-field strength into a lower 
sheet, where (B 2 ) 112 increases as a function of the magnetisation 
parameter Aq, and an upper sheet, where the opposite happens. The 
left panel of Fig. [9] in particular, shows the surface deformation 
e s for the lower part of the solution space, while the right panel is 
its continuation beyond the turning point of (B 2 ) 1 ' 2 and thus repre- 
sents the upper sheet of the surface. Since nonrotating magnetised 
models always have a prolate shape, i.e. e s < 0, while rotating un- 
magnetised models always have an oblate one, i.e. e s > 0, it follows 
that rotating models between these limiting cases are divided into 
prolate ones and oblate ones by a neutral line e s = 0. According to 
its definition in Eq. l |68[ >, the neutral line e s = does not require the 
stellar interior to be spherically symmetric and, indeed, when mov- 
ing towards large values of the circumferential radius B circ , models 
with comparable values of the equatorial and polar coordinate radii 
r e , r p (i.e. with e s — > 0), show a progressively more pronounced 
"diamond shape" caused by an increasingly spindle-shaped matter 
distribution inside the star. 

Unlike the magnetic potential, the centrifugal one is not con- 
fined to the star, and its influence increases moving away from the 
axis of rotation, reaching its minimum at the equator of the star. 
Note that M = at the surface of the star, so that the latter is now 
an isosurface of the function v - lnT [cf. Eq. l|53[>]. From Fig. [9] 
we further infer that for any fixed angular velocity SI > Sl , where 
flo = 2.594 x 10 3 s _1 is the maximum angular velocity at the mag- 
netisation truncation limit, the mass-shedding limit is reached for 
some oblate shape and e s > 0. Since the toroidal magnetic field 
is a source of additional pressure, augmenting the magnetisation 
not only increases the deformation of the surface and of the mat- 
ter distribution, but it also causes an expansion, particularly of the 
circumferential radius B C j rc . At a sufficient level of magnetisation, 
it will therefore be possible to reach the mass-shedding angular ve- 
locity, i.e. the rotation frequency such that the condition of geodesic 
motion at the equator is satisfied, and the star will develop the 
characteristic cusp at the equator. Because this angular frequency 
is smaller than the corresponding one for an unmagnetised model 
having the same rest-mass, the toroidal magnetic field indirectly 
sets a reduced limit of the spin frequency of these objects. In Fig. [9] 
this mass-shedding limit corresponds to the lower right boundary of 
the upper part of the solution space (see right panel) and along this 
limit, the Newtonian condition of geodesic motion, QrR 3 = const, 
is fulfilled. Note also that with the exception of the unmagnetised 
one, all mass-shedding models appear to belong to the upper sheet 
of the solution space. For models rotating at angular velocities up 
to that of model 00-Pol2, i.e. SI = 5.050 x 10 3 s~\ we have ver- 
ified this proposition directly by determining the turning point of 
(B 2 } 1/2 ; for even more rapidly rotating models for which the nu- 
merical determination of the maximum field-strength limit has not 
been conclusive, it is supported by extrapolating the boundary be- 
tween lower and upper part of the solution space beyond model 
00-Pol2 towards the unmagnetised mass-shedding limit. 



Using a small set of additional rotating models computed be- 
yond the magnetisation truncation limit, we have started to explore 
the behaviour of the equilibrium models in these rather extreme 
conditions. Overall, we have found that the angular velocity of 
mass-shedding models decreases progressively, while the bound- 
ary associated with the mass-shedding limit and the neutral line 
6 S = converge to the point (SI 2 , (B 2 )) = (0, 0). These results sug- 
gest therefore that the solution at (Si 2 , (B 2 )) = (0, 0) in the upper 
sheet of the space of solutions corresponds to the limit of a nonro- 
tating model of infinite radius and vanishing mean magnetic-field 
strength. This fascinating suggestion clearly requires a more ex- 
tensive analysis to be confirmed; we postpone this to a subsequent 
work. 

Figure [5] also reveals that unlike in the unmagnetised rotat- 
ing case with the same angular velocity SI, all magnetised mass- 
shedding configurations behave qualitatively identically when the 
magnetisation is altered, namely: lowering their magnetisation 
moves them away from the mass-shedding limit. Accordingly, rota- 
tion becomes sub-critical, and the characteristic cusp at the equator 
disappears, which is accompanied by a decrease of the surface de- 
formation. Slowly-rotating models pass through different stages as 
the magnetisation is lowered from the mass-shedding limit to the 
given angular velocity SI down to the limit of vanishing magneti- 
sation. More specifically, we find that: (i) below the mass-shedding 
limit, 6 S remains positive, but decreases continuously until the neu- 
tral line e s = is crossed for the first time; (ii) models then be- 
come increasingly prolate because the toroidal magnetic field be- 
comes the principal source of deformation until they reach some 
negative minimum value of e s ; (iii) as the magnetisation is fur- 
ther decreased, e s increases progressively, and models eventually 
become oblate again when they cross the neutral line e s = for a 
second time, and rotation prevails over a decreasing toroidal mag- 
netic field. For models rotating at moderate angular velocities of 
3 x 10 3 s~' < SI < 5 x 10 3 s -1 , the behaviour is similar to that 
of slowly rotating ones, except that e s remains positive, so that the 
models always exhibit an oblate shape regardless of the level of 
magnetisation. Finally, rapidly rotating models with SI > 5x 10 3 s _I 
remain close to the magnetised mass-shedding limit regardless of 
the magnetisation level. 

Generally speaking, for any fixed angular velocity SI, the line 
SI = const is tangential to exactly one level curve of e s . The touch- 
ing point determines the magnetisation level for which the surface 
deformation attains its minimum. For decreasing angular velocity, 
these minima move to higher magnetisations and smaller values of 
e s . The angular velocity SI' for which the maximum magnetic-field 
strength model exhibits e s = separates rotating models into two 
groups: (1) models rotating at SI > SI' are always oblate, and the 
model of minimum (positive) surface deformation is located in the 
lower part of the solution space; (2) models rotating at SI < SI' in- 
clude a prolate model of minimum (negative) surface deformation 
which is located in the upper part of the solution space. 

Figure [TO] provides an equivalent representation of the data 
shown in Fig. [9] but this time in terms of the quadrupole distor- 
tion e. Also in this case, the space of the numerical solutions is 
split into an upper sheet and a lower one, where the distinction is 
the same one as made for the surface distortion e s . For all models, 
6 increases monotonically as a function of the magnetisation, and 
the neutral line e = extends from the nonrotating and unmagne- 
tised reference model, up to a strongly-magnetised mass-shedding 
model with (B 2 ) 1 ' 2 s 2x 10 17 G rotating at an angular velocity 
of SI « 4.5 x 10 3 s -1 , and with T/\W\ « 0.06. Thus, even rapidly 
rotating models can exhibit a prolate matter distribution provided 
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Table 4. Neutron star models at the maximum field-strength limit. The gravitational mass is M = 1 .400 Mq for each EOS respectively in the unmagnetised 
and nonrotating case for which properties are listed in Table [5] B max is the maximum value of the magnetic-field strength, (B 2 ) 1 ' 2 the root mean square 
value of the magnetic-field strength determined over the volume of the star, ^/\W\ the absolute value of the ratio of total magnetic energy ^# and potential 
energy W, «b c the central baryon number density [0.1 fin ], M the gravitational mass, i? c j lc the circumferential radius, / the moment of inertia, e s the surface 
deformation, eg the magnetic quadrupole distortion, and GRV2/GRV3 the estimates of the global error of respective models based on the relativistic virial 
identities introduced in Sect. 14. 31 
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the magnetisation is sufficiently strong. Although we have trun- 
cated the solution space toward high magnetisations, the minimum 
quadrupole distortion of e = -0.6127, which corresponds to the 
strongest prolate deformation of our sample and has been obtained 
for the nonrotating model at the magnetisation truncation limit, ex- 
ceeds considerably the quadrupole distortion e = 0.1537 of the un- 
magnetised mass-shedding model rotating at Q = 5.205 x 10 3 s 
with T/\ W\ = 0.09380. This significant difference can be explained 
by the large ratio = 0.2448 of the magnetised model. Since 

no mass-shedding limit has been encountered in the nonrotating 
magnetised case, ^ l\W\ can possibly attain arbitrarily high values 
and accordingly, the quadrupole distortion e can, in principle, grow 
without bounds, too. We leave the assessment of this conjecture to 
a subsequent work. 

We also note that whereas all mass-shedding models exhibit 
an oblate surface deformation for which e s > 0, only rapidly rotat- 
ing mass-shedding models with £2 > 4.5 x 10 3 s show a positive 
quadrupole distortion with e > 0. On the other hand, all mass- 
shedding models rotating at lower angular velocities actually pos- 
sess a prolate matter distribution and e < which, moreover, seems 
to grow without bounds for increasing magnetisation like in the 
nonrotating case. Altogether, Figs.|9]and|10|reveal that the neutral 
lines 6 S = and e = differ significantly, suggesting the division 
of the space of solutions of magnetised and rotating models into 
three classes: (1) models for which apparent shape and distortion 
of the matter distribution are both prolate, thus with 6 S < and 
6 < 0, which we label PP for prolate-prolate; (2) models for which 
the apparent shape is oblate and the distortion of the matter distri- 
bution is prolate, thus with e s > and e < 0, which we label PO 
for prolate-oblate; (3) models for which apparent shape and distor- 
tion of the matter distribution are both oblate, thus with e s > and 
e > 0, which we label 00 for oblate-oblate. As a result, in contrast 
with the results of Kiuchi & Yoshida ( 2008 1, rotating models with a 
strong toroidal magnetic field do not necessarily exhibit a negative 
quadrupole distortion e [cf. discussion in Sect.|2]and Fig.^. 

Representative models for each of the three classes located at 
the maximum field-strength limit have been presented in Sect. [5] 
namely model PP-Pol2, the nonrotating configuration with a pro- 
late apparent shape and a prolate matter distribution, and earlier 
in this Section, namely model P0-Pol2, with a prolate apparent 
shape and an oblate matter distribution, and model 00-Pol2 with 
an oblate shape and an oblate matter distribution. Within Figs. [5] 
and [TO] their positions in the lower and upper parts of the solution 
space have been marked by red dots along the line of maximum 



magnetic-field strength related to (B 2 ) 1 ' 2 . Note that no model ex- 
ists for which the apparent shape is prolate and the distortion of the 
matter distribution oblate, i.e. e s < and e > 0; therefore, because 
of the different nature of the forces caused by a toroidal magnetic 
field and rotation, a class OP does not exist. While magnetic and 
centrifugal forces distort the matter distribution at a comparable 
level with respect to the involved amounts of magnetic and kinetic 
energy, the magnetic potential is confined and does not act on the 
surface of the star directly through Eq. j53}. Therefore, already in 
the magnetised and nonrotating case, the surface deformation is al- 
ways smaller than the quadrupole distortion. Moreover, since cen- 
trifugal forces act more efficiently at large distances from the rota- 
tion axis, their influence on the surface of the star can be significant, 
and indeed the value of e s = 0.7242 obtained at the unmagnetised 
mass-shedding limit is considerably larger (in absolute terms) than 
the value obtained for the nonrotating model at the upper magneti- 
sation limit of e s = -0.1369; the opposite is true when considering 
the quadrupole distortion with respective values of e = -0.6127 
and 6 = 0.1537. 

In order to assess the general validity of our results obtained 
in the rotating case for the analytic Pol2 EOS, we have compared 
the two rotating models P0-Pol2 and 00-Pol2 with the corre- 
sponding models PO-APR and 00-APR, built with the APR EOS 
for identical values of the T/\W\ ratio at the unmagnetised limit, 
i.e. T/\W\ = 0.03 and T/\W\ = 0.07. The physical properties of all 
models are listed in Table [4] where it is clear that the effects of ro- 
tation are less pronounced for the two APR models because of the 
stiff er nature of the APR EOS. The APR models admit a higher 
maximum value of T/\W\ = 0.1052 at the unmagnetised mass- 
shedding limit, which should be compared with the corresponding 
T/\W\ = 0.09380 for the softer Pol2 EOS. Moreover, as discussed 
in Sect. [5] the APR EOS reference model supports a higher max- 
imum field strength when compared to the Pol2 EOS one, which 
causes magnetic effects to be stronger for the two models PO-APR 
and 00-APR. 

The toroidal magnetic-field strength and the baryon number 
density for the PO-APR and 00-APR models are reported in Fig. 1 1 1 1 
which shows that the matter distributions are less condensed (the 
EOS is stiffer) and the peaks of the magnetic-field strength have 
moved slightly outward. Their prolate deformation appears more 
pronounced in agreement with the higher ratio *j?/\W\, and the 
outer crust is easily discernible because the magnetic field is es- 
sentially absent from this low density region. The ratio T/\W\ for 
models PO-APR and 00-APR is smaller than that of their unmag- 
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Figure 11. Isocontours of magnetic-field strength (left panel) and baryon number density (right panel) in the (x, z) plane of the PO-APR and 00-APR models of a 
star built with the APR EOS with a gravitational mass of M = 1.400Mq in the unmagnetised and nonrotating case which is now rotating at CI = 4.2 19 X 10 3 s _1 
(top) and fi = 6.004X 10 3 s (bottom) respectively. The gray disc indicates the dimensions of the unmagnetised reference model. Physical properties of models 
PO-APR and 00-APR are listed in Tableg] 



Table 5. Static reference models without magnetic field and with a gravitational mass of M = 1 .400 Mq . nt, rQ is the central baryon number density [0. 1 fm ] , 
Mb the baryon mass, R aic the circumferential radius, / the moment of inertia, bg the magnetic distortion coefficient defined in Eq. (73) , the rotational 
distortion coefficient defined in Eq. (73) , eg the magnetic distortion coefficient defined in Eq. (73) , cq the rotational distortion coefficient defined in Eq. (73) , 
and GRV2/GRV3 the estimates of the global error of respective models based on the relativistic virial identities introduced in Sect.|4.3| 
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Figure 12. Surface deformation e s (left panel) and quadrupole distortion e (right panel) for rotating models built with the Pol2 EOS and the APR EOS as a 
function of the mean square magnetic-field strength (B 2 ) from the unmagnetised limit up to the maximum field strength models P0-Pol2 and PO-APR as well 
as 00-Pol2 and 00-APR. The angular velocities are chosen such that T/\W\ = 0.03 and T/\W\ = 0.07 respectively in the unmagnetised case. 



netised counterparts, in agreement with their smaller momenta of 
inertia. On the contrary, models P0-Pol2 and 00-Pol2 are already 
located in the region of increasing momenta of inertia and show 
ratios TJ\W\ = 0.03200 and T/\W\ = 0.07119, which are larger 
than those of the unmagnetised models rotating at the same angular 
velocities SI. 

Finally, Fig. [12] compares the dependence of e s and e on 
the magnetisation levels for the four reference rotating models 
(cf. Fig. [5] where the comparison was made for the nonrotating 
models PP-Pol2 and PP-APR). Clearly, there is a good qualita- 
tive agreement between the two EOSs with differences that are 
due mostly to the higher maximum magnetic-field strength and the 
higher maximum value of T/\W\ supported by models PO-APR and 
00-APR. Note that model 00-Pol2 is so close to the mass shedding- 
limit that it is the only one for which models rotating at the same 
angular velocity SI show increasing e s from the unmagnetised limit 
on. While all sequences in the left panel of Fig. [T2] maintain an 
oblate shape, those associated with models P0-Pol2 and PO-APR 
show a transition from an oblate matter distribution to a prolate one 
(right panel of Fig.[T2|. 



7 DISTORTION COEFFICIENTS 

Despite the complex behaviour shown by the equilibrium models 
when both the magnetic-field strength and the rotation are varied, it 
is possible to express such a behaviour through a very simple alge- 
braic expression. This was pointed out already by Wentzel ( 1960) 
and by |Ostriker & GurrR] ( |1969] l, who have considered this issue 
in earlier Newtonian studies and have suggested to parametrize the 
quadrupole distortion e Newt induced by a toroidal magnetic field and 
by rotation in a self-gravitating incompressible fluid respectively as 



JK T 
e N e„. - cb + en - -a B — + a a — 



(72) 



where a B = a n = 3.750. This approximation was adopted also 
by Cutler ( 2002 ) in order to derive an estimate for the quadrupole 
distortion of neutron stars, within a Newtonian framework. Since 
neutron stars are highly relativistic objects, we next consider 
whether an expression similar to Eq. \72\ can be derived in a rel- 
ativistic regime and the quantitative differences that then emerge 
with respect to a Newtonian treatment. 



We start by recalling that already in Sect. [6] we have remarked 
how the surface deformation e s and the quadrupole distortion e for 
highly magnetised and rapidly rotating models exhibit an almost 
linear dependence on SI 2 and on (B 2 ). As a result, we can express 
these quantities in terms of "deformation coefficients" 



6 S = -b B {B 2 l5 ) + b n Sl 2 



i <flj 5 > + c a si 2 



(73) 



where S 15 expresses the magnetic-field strength in units of 10 15 G 
and SI is expressed in s _I . The distortion coefficients b B , b& for the 
surface deformation e s and the coefficients c B , c n for the quadrupole 
distortion e can be easily computed through the directional deriva- 
tives along the coordinate axes of SI 2 and {B 2 ) and are collected 
in Table [5] besides the basic properties of the unperturbed mod- 
els for all of the EOSs considered. The importance of these dis- 
tortion coefficients is that they provide all the information needed 
to compute the resulting values of e s and of e by simply insert- 
ing appropriate values for the angular velocity SI and for the mean 



magnetic-field strength (Bf 5 ) ' in Eq. 1 73 1. In addition, the approx 



imate expressions of Eq. {73} can cover a very large portion of the 
physically-realistic space of parameters and, for example, in the 
case of the Pol2 EOS the phenomenological relations yield relative 
errors < 8% up to values of (B 2 ) [/2 = 1 x 10 17 G and angular veloc- 
ities of SI = 2 x 10 3 s , which largely cover all known magnetars. 

We note that the quadratic dependence expressed by Eq. l |73| > 
applies only to old and cold neutron stars without a superconduct- 
ing proton phase, where the latter is expected to be suppressed 
for very strong magnetic fields above the critical one, i.e. for 
(5) > fi c i s lx 10' 5 G. On the other hand, for magnetic field 
strengths of (B) < B cl , superconductivity is possible and the mag- 
netic field is confined to flux tubes of strength 6 C [ « 1 x 10 13 G, 
which enlarge the anisotropic part of the average electromagnetic 
stress tensor (Jones 1975, Easson & Pethick 1977), leading to a 
modification of the magnetic forces. As a result, for magnetised 
models with {B} < 10' 5 G, a linear dependence of the matter distor- 
tion on the average magnetic-field strength is expected and Eq. l |73| > 
has to be modified taking the alternate form 



e s = -b B (B ls ) + b n Sl 2 



:(B l5 ) + c a SI 2 



(74) 



where the coefficients in Eq. l |74[ > are the same as those in Eq. (73} . 

Figure [T3] offers a graphical representation of the distortion 
coefficients for the surface deformation (left panel) and for the 
quadrupole distortion (right panel). The coefficients b B and c B 
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Figure 13. Distortion coefficients fog and b& for e s (left panel) as well as eg and cq for e (right panel) derived in the linear regime by perturbing a nonrotating 
and unmagnetised star with a gravitational mass of M = 1.400 Mq. 



as well as b a and c a are computed for different EOSs and us- 
ing the reference nonrotating models with a gravitational mass of 
M = 1.400 Mq. Also shown as a gray-shaded band are the dis- 
tortion coefficients of nonrotating models with a Pol2 EOS hav- 
ing the same gravitational mass, but where the polytropic constant 
is varied to obtain different circumferential radii (these have 
therefore R cilc 3> 9.891 km). The influence of the different EOSs is 
in fact most easily followed through the associated circumferential 
radius R CKC , which had already turned out to be a crucial quantity 
for the maximum field strength models presented in Sect. [5] As 
we explain in Appendix|A] using basic scaling considerations, the 
magnetic distortion coefficients exhibit a dependence of the type 
bn.Cn oc R 4 . , whereas the rotational distortion coefficients exhibit 

a » a clrc s 

a dependence of the type bn, cq oc R\ 

Also shown for comparison in Fig.[T3]are the Newtonian ref- 
erence model from Cutler (2002), labelled Pol2C18, and its numer- 
ical equivalent as computed with our code in the Newtonian limit, 
labelled as Pol2N160 Note that both Newtonian models have sim- 
ilar values of c B and that the Pol2N18 data point is rather close 
to the reference band of relativistic Pol2 EOS models. However, 
the rotational distortion coefficient c a reported in Cutler ( 2002) ex- 
ceeds the correct Newtonian one of model Pol2NlQ by almost a 
factor of two. Furthermore, both data points largely overestimate 
the correct relativistic result obtained for model Pol2R10 which 
is about a factor of four smaller and located at the lower end of 
the relativistic Pol2 EOS. This simple example highlights therefore 
how a Newtonian treatment is inadequate for the determination of 
realistic distortion coefficients for real neutron stars. 



8 CONCLUSIONS 

We have computed models of rotating relativistic stars with a 
toroidal magnetic field under the assumption that the matter is a 
single-constituent perfect fluid described by a one-parameter equa- 
tion of state and behaves as a perfect conductor subject to the 
laws of ideal magnetohydrodynamics. We have investigated the 
combined effects of a toroidal magnetic field and rotation on the 
apparent shape and on the internal matter distribution, focussing 
in particular on the quadrupole distortion, as this is the relevant 



6 Only the quadrupole distortion coefficients eg and cq for model Pol2C18 
were considered by Cutler 1 2002j. 



quantity for the gravitational-wave emission. Models of maximum 
field strength have been computed for a sample of eight different 
nuclear-matter equations of state, together with the surface defor- 
mation and the quadrupole distortion. 

We have found that nonrotating models appear to admit arbi- 
trary levels of magnetisation accompanied by a seemingly unlim- 
ited growth of size and quadrupole distortion. In particular, we have 
been able to compute a highly magnetised model for the Pol2 EOS 
with a baryon mass of M b = 1.680 Mq, whose circumferential ra- 
dius of /? c ; rc = 14.30 km in the unmagnetised case, inflates to i? C i rc = 
101.5 km for an average magnetic field of (B) = 0.1461 x 10 17 G. 
These results should be contrasted with those of Kiuchi & Yoshida 
( 2008), who reported a loss of convergence for this model at a mod- 
erate level of magnetisation and corresponding to a value of merely 
/f circ = 28.85 km. 

When considering rotating models we have instead found that 
the increase in equatorial size introduced by the toroidal magnetic 
field reduces the frequency at which mass-shedding would other- 
wise appear in unmagnetised models. Overall, the full space of so- 
lutions can be split up into three distinct classes for which the sur- 
face distortion and the quadrupole distortion are either prolate and 
prolate, oblate and prolate or oblate and oblate, respectively. 

We have also determined the relativistic distortion coefficients 
whose absolute value depends mainly on the radius of the star for 
all the EOSs considered. Using such coefficients it is possible to 
compute the surface deformation and the quadrupole distortion by 
means of a simple algebraic expression which is effective for all 
magnetisations and rotation rates of known magnetars. Finally, a 
comparison with the corresponding Newtonian distortion coeffi- 
cients has shown that the latter overestimate the quadrupole dis- 
tortion induced by the toroidal magnetic field by about a factor 
of five and the one induced by rotation by about a factor of four. 
Hence, they are inadequate for strongly-relativistic objects like neu- 
tron stars. 

The results presented here relative to equilibrium configura- 
tions provide the first basic steps to explore the stability properties 
of magnetised stars, whose analysis in full general relativity has re- 
cently seen a spur of activity ( [Kiuchi et al.|201 1| Lasky et al.|201 1| 
|Ciolfi et al-lpOTT] [ZrnFet al.||2012| |Lasky et al.||2012| |Ciolfi &] 
Rezzolla 2012). We will investigate the stability properties of our 
models in a forthcoming future work. 
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APPENDIX A: ASSESSMENT OF RESULTS BY CUTLER 
(2002) 

As anticipated in Sect. [7] earlier Newtonian studies have sug- 
gested to express the quadrupole distortion e Ncwt induced in a self- 
gravitating incompressible fluid by a toroidal magnetic field and by 
rotation respectively as function of the total magnetic energy ', of 
the kinetic energy T, and of the potential energy W, namely through 
Eq. {72) < |Wentzel|1960||Ostriker & Gunn|[T 969). For a spherical 
star built with a polytropic EOS with y = 2 it is easy to estimate 
that 

T=-k 1 MR 2 £1 2 , W = --—, (Al) 
5 4 R 

where K[ = 0.65345 is a constant derived by |Lai et al.|(1993) . As 
a result, the ratio of the kinetic to binding energy will scale as 
T l\W\ oc R 3 . On the other hand, the total magnetic energy ^ of 
the same body reads 

i r i 47: , , 

JC= — B 2 dV = —R 3 (B 2 ), (A2) 
Sn J v 3 

where (B 2 ) denotes the mean square average of the magnetic-field 
strength B inside the star. It then follows that ^/\W\ oc R 4 . Even- 
tually, the total distortion e can be expressed in terms of mean 
magnetic-field strength (B 2 ) 1 ' 2 and angular velocity £2, namely 
[cf.Eq.j73j] 

fNcwt = e b + £r = -c B (B 2 15 ) + c a Q 2 . (A3) 

For a Newtonian model with M = 1 .400 Mq and a radius of R = 
10.00km, Cutler (2002j has computed c B according to Eqs. \Al\ 
and {A2} and reported the distortion coefficients as 

c„ = 1.600 x 10~ 6 , c a = 7.600 x 10~ 9 . (A4) 

However, if we use Eqs. ( |Alfr and \A2\ and adopt an identical 
model with M = 1 .400 Mq and R = 10.00 km, we obtain 

c a = 1.610x10-*, c n = 3.516 x 10~ 9 , (A5) 

where c n = 6.725 x 10~ 9 for an incompressible fluid. Although 
the estimates for c B agree reasonably well, our value for c n is less 
than half the one quoted in |Cutler| ( |2002| >, which is instead close to 
the corresponding value for a homogeneous sphere. To validate our 
estimates, we have additionally computed the distortion coefficients 
with our numerical code in the Newtonian limit and obtained 

c B = 1.51 1 x 10- 6 , c a = 3.569 x 10~ 9 , (A6) 

which agree well with the estimates from Eq. l|A5j. The remain- 
ing difference is due to the fact that a B = a a = 3.750 only in the 
incompressible case and need to be corrected in the compressible 
one. After taking into account this EOS dependence, full agreement 
is achieved, as can be verified in Table |Al| which provides a com- 
pilation of the different values for the coefficients c B and c n and 
where the coefficients a B and an have been added when available. 



Table Al. Distortion coefficients for a Newtonian star with M = 1 .400 Mq 
and R = 10.00 km for a polytropic EOS with y = 2. For comparison, values 
for an incompressible model with the same properties are included. 



r 


as 


aa 


CB 




ca 




oo 1 


3.750 


3.750 


2.013 x 10' 


-6 


6.725 x 10" 


-9 


2 2 


3.750 


3.750 


1.610 x 10' 


-6 


3.516 x 10" 


-9 


V 


3.518 


3.804 


1.511 x 10" 


-6 


3.569 x 10" 


-9 


2 4 






1.511 x 10" 


-6 


3.567 x 10" 


-9 


2 5 


3.750 




1.600 x 10" 


-6 


7.600 x 10" 


-9 



1. Distortion coefficients c B and cq according to Eq. j72) . In the incom- 
pressible case, T = I£l 2 /2 with / = (2/5)MR 2 , and W = -(3/5)M 2 /«. 

2. As before but using Eqs. {AT} and (A2) . 

3. Distortion coefficients derived from the linear perturbation of a numerical 
Newtonian model. 

4. Distortion coefficients from Haskell et al. 1 2008) revised by Hask ell et al.| 
(2009) . 

5. Distortion coefficients from Cutler 1 2002 1. 
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